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O : ABSTRACT 

b 

In this paper, we consider minimizing the action functional as a method for 
numerically discovering periodic solutions to the n-body problem. With this 
method, we can find a large number of choreographies and other more general 
solutions. We show that most of the solutions found, including all but one of 
the choreographies, are unstable. It appears to be much easier to find unstable 
solutions to the n-body problem than stable ones. Simpler solutions are more 

CT} ■ likely to be stable than exotic ones. 
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1. Least Action Principle 



2 ■ Given n bodies, let rrij denote the mass and Zj(t) denote the position in R 2 = C of 

body j at time t. The action functional is a mapping from the space of all trajectories, 
zi(t), z 2 (t), . . . , z n (t), < t < 2tt, into the reals. It is defined as the integral over one period 
of the kinetic minus the potential energy: 
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Stationary points of the action function are trajectories that satisfy the equations of 
motions, i.e., Newton's law gravity. To see this, we compute the first variation of the action 
functional, 
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and set it to zero. We get that 



zf-z: 



" E m J m ^ lis' .7 = !,2,...,n, a = 1,2 

^— ' I ! ■ — - , 1 1 •> 



I % Zk | 



Note that if rrij = for some j, then the first order optimality condition reduces to 
= 0, which is not the equation of motion for a massless body. Hence, we must assume that 
all bodies have strictly positive mass. 



2. Periodic Solutions 

Our goal is to use numerical optimization to minimize the action functional and thereby 
find periodic solutions to the n-body problem. Since we are interested only in periodic 
solutions, we express all trajectories in terms of their Fourier series: 

oo 
k=— oo 

Abandoning the efficiency of complex-variable notation, we can write the trajectories with 
components Zj(t) = (xj(t),yj(t)) and ^ k = (ak,Pk)- So doing, we get 

00 

x(t) = <2q + (a c k cos(kt) + a s k sm(kt)) 

fc=i 

00 

y(t) = b + J2( b tcos(kt)+b s k sm(kt)) 
k=i 

where 

Oo — Oio, a c k = a>k + oi-k-, of k = [3_k ~ Pk, 
b = fa, b c k = (3 k + P-k, K = a k~ a-fc- 

Since we plan to optimize over the space of trajectories, the parameters do, a c k , a s k , bo, b k , 
and b k are the decision variables in our optimization model. The objective is to minimize 
the action functional. 

ampl is a small programming language designed for the efficient expression of optimiza- 
tion problems ?. Figure 1 shows the AMPL program for minimizing the action functional. 

Note that the action functional is a nonconvex nonlinear functional. Hence, it is expected 
to have many local extrema and saddle points. We use the author's local optimization 
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software called LOQO (see ?, ?) to find local minima in a neighborhood of an arbitrary 
given starting trajectory. One can provide either specific initial trajectories or one can give 
random initial trajectories. The four lines just before the call to solve in Figure 1 show 
how to specify a random initial trajectory. Of course, AMPL provides capabilities of printing 
answers in any format either on the standard output device or to a file. For the sake of 
brevity and clarity, the print statements are not shown in Figure 1. AMPL also provides the 
capability to loop over sections of code. This is also not shown but the program we used 
has a loop around the four initialization statements, the call to solve the problem, and the 
associated print statements. In this way, the program can be run once to solve for a large 
number of periodic solutions. 

2.1. Choreographies 

Recently, ? introduced a new family of solutions to the n-body problem called chore- 
ographies. A choreography is defined as a solution to the n-body problem in which all of 
the bodies share a common orbit and are uniformly spread out around this orbit. Such 
trajectories are even easier to find using the action principle. Rather than having a Fourier 
series for each orbit, it is only necessary to have one master Fourier series and to write the 
action functional in terms of it. Figure 2 shows the AMPL model for finding choreographies. 

3. Stable vs. Unstable Solutions 

Figure 3 shows some simple choreographies found by minimizing the action functional 
using the AMPL model in Figure 2. The famous 3-body figure eight, first discoverd by ? 
and later analyzed by ?, is the first one shown — labeled FigureEight3. It is easy to find 
choreographies of arbitrary complexity. In fact, it is not hard to rediscover most of the 
choreographies given in ?, and more, simply by putting a loop in the AMPL model and 
finding various local minima by using different starting points. 

However, as we discuss in a later section, simulation makes it apparent that, with the 
sole exception of FigureEight3, all of the choreographies we found are unstable. And, the 
more intricate the choreography, the more unstable it is. Since the only choreographies that 
have a chance to occur in the real world are stable ones, many cpu hours were devoted to 
searching for other stable choreographies. So far, none have been found. The choreographies 
shown in Figure 3 represent the ones closest to being stable. 

Given the difficulty of finding stable choreographies, it seems interesting to search for 
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stable nonchoreographic solutions using, for example, the ampl model from Figure 1. The 
most interesting such solutions are shown in Figure 4. The one labeled Ducati3 is stable 
as are Hill3_15 and the three DoubleDouble solutions. However, the more exotic solutions 
(OrthQuasiEllipse4, Rosette4, PlateSaucer4, and BorderCollie4) are all unstable. 

For the interested reader, a JAVA applet can be found at ? that allows one to watch the 
dynamics of each of the systems presented in this paper (and others). This applet actually 
integrates the equations of motion. If the orbit is unstable it becomes very obvious as the 
bodies deviate from their predicted paths. 



3.1. Ducati3 and its Relatives 

The Ducati3 orbit first appeared in ? and has been independently rediscovered by this 
author, Broucke ?, and perhaps others. Simulation reveals it to be a stable system. The 
JAVA applet at ? allows one to rotate the reference frame as desired. By setting the rotation 
to counter the outer body in Ducati3, one discovers that the other two bodies are orbiting 
each other in nearly circular orbits. In other words, the first body in Ducati3 is executing 
approximately a circular orbit, zi(t) = —e lt , the second body is oscillating back and forth 
roughly along the rr-axis, ^(t) = cos(t), and the third body is oscillating up and down the 
?/-axis, z 3 (t) = isin(t). Rotating so as to fix the first body means multiplying by e~ lt : 

Zl (t) = e - u (-e u ) = -l 

z 2 (t) = e~ lt cos(t) = (l+e~ 2it )/2 

z 2 (t) = e- u ism(t) = (1 -e- 2it )/2. 

Now it is clear that bodies 2 and 3 are orbiting each other at half the distance of body 1. 
So, this system can be described as a Sun, Earth, Moon system in which all three bodies 
have equal mass and in which one (sidereal) month equals one year. The synodic month is 
shorter — half a year. 

This analysis of Ducati3 suggests looking for other stable solutions of the same type but 
with different resonances between the length of a month and a year. Hill3_15 is one of many 
such examples we found. In Hill3_15, there are 15 sidereal months per year. Let Hill3_n 
denote the system in which there are n months in a year. All of these orbits are easy to 
calculate and they all appear to be stable. This success suggests going in the other direction. 
Let Hill3_^- denote the system in which there are n years per month. We computed Hill3_| 
and found it to be unstable. It is shown in Figure 6. 

In the preceding discussion, we decomposed these Hill-type systems into two 2-body 
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problems: the Earth and Moon orbit each other while their center of mass orbits the Sun. 
This suggests that we can find stable orbits for the 4-body problem by splitting the Sun into 
a binary star. This works. The orbits labeled DoubleDoublen are of this type. As already 
mentioned, these orbits are stable. 

Given the existence and stability of FigureEight3, one often is asked if there is any 
chance to observe such a system among the stars. The answer is that it is very unlikely 
since its existence depends crucially on the masses being equal. The Ducati and Hill type 
orbits, however, are not constrained to have their masses be equal. Figure 5 shows several 
Ducati-type orbits in which the masses are not all equal. All of these orbits are stable. This 
suggests that stability is common for Ducati and Hill type orbits. Perhaps such orbits can 
be observed. 



4. Limitations of the Model 

The are certain limitations to the approach articulated above. First, the Fourier series 
is an infinite sum that gets truncated to a finite sum in the computer model. Hence, the 
trajectory space from which solutions are found is finite dimensional. 

Second, the integration is replaced with a Riemann sum. If the discretization is too 
coarse, the solution found might not correspond to a real solution to the n-body problem. 
The only way to be sure is to run a simulator. 

Third, as mentioned before, all masses must be positive. If there is a zero mass, then 
the stationary points for the action function, which satisfy (1), don't necessarily satisfy the 
equations of motion given by Newton's law. 

Lastly, the model, as given in Figure 1, can't solve 2-body problems with eccentricity. 
We address this issue in the next section. 



5. Elliptic Solutions 

An ellipse with semimajor axis a, semiminor axis b, and having its left focus at the 
origin of the coordinate system is given parametrically by: 

x(t) = f + a cost, y(t)=bsint, 

where / = \/a? — b 2 is the distance from the focus to the center of the ellipse. 

However, this is not the trajectory of a mass in the 2-body problem. Such a mass 
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will travel faster around one focus than around the other. To accomodate this, we need to 
introduce a time-change function 9{t): 

x(t) = f + a cos 9 (t), y(t) = bsm9(t). 

This function 9 must be increasing and must satisfy 9(0) = and 9(2n) = 2ir. 

The optimization model can be used to find (a discretization of) 9(t) automatically by 
changing param theta to var theta and adding appropriate monotonicity and boundary 
constraints. In this manner, more realistic orbits can be found that could be useful in real 
space missions. 

In particular, using an eccentricity e = f /a = 0.0167 and appropriate Sun and Earth 
masses, we can find a periodic Hill- Type satellite trajectory in which the satellite orbits the 
Earth once per year. 

6. Sensitivity Analysis 

The determination of stability vs. instability mentioned the previous sections was done 
empirically by simulating the orbits with a integrator and very small step sizes. Two inte- 
grators were used: a midpoint integrator and a 4-th order Runge-Kutta integrator. Orbits 
that are claimed to be stable were run for several hours of cpu time (which corresponds to 
many thousands of orbits) without falling apart. Orbits that are claimed to be unstable 
generally became obviously so in just a few seconds of cpu time, which corresponds to only 
a few full orbits. In this section, we describe a Floquet analysis of stability and present this 
measure of stability for the various orbits found. 

For simplicity, in this section we assume that all masses are equal to one. Let £*(t) = 
(z*(t), z*(t)) be a particular solution to 

i = a(o 

where 

A(z(t),z(t)) = (z(t),a(z(t))) 

and 

a(z) = (ai(z), . . .,a n (z)) 

and 

a A z ) = ~ J = l,2,...,n. 
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Consider a nearby solution £(£): 



i(t) = Am) 

« A(r(t)) + A'(rw)(ew-rw) 
= r(t) + A'(r(t))(ew-rw). 

Put A£ = £ - £*. Tnen A£ = A finite difference approximation yields 

A£(t + h) = A£(t) + hA\Z*(t))AZ(t) 
= (I + hA'(C(t)))m)- 

Iterating around one period, we get: 



'n-\ 



where h — T/n and U = iT/n. 

The following perturbations, which are associated with invariants of the physical laws, 
are unimportant in calculating A£(T): 
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where i? denotes rotation by 90°. The first two of these perturbations correspond to trans- 
lation. The next two correspond to moving frame of reference and the last two correspond 
to rotation, and dilation. Dilation is explained below. Of course, all positions and velocities 
are evaluated at t — 0. Vector d denotes the i-th unit vector in M 2 . 

Consider spatial dilation by p together with a temporal dilation by 9: 

Z J (t)=pz J (t/6). 
Given that the solution, it is easy to check that 

Zj(t)-Z k (t) 



Mi 



' \\Zj(t) - Z k (tW 
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Hence, if mass is to remain fixed, we must have that p 3 = 9 2 : 

Zj(t) = pZj (t/p 3/2 ) z s {t) = p-WzMf?' 2 ). 

To find the perturbation direction corresponding to this dilation, we differentiate with respect 
to p at p = 1 : 



d 

dp 



pzjit/p^) 

Zj(t/ff 



-1/2 i 



3/2A 



p=l 



' 2 z i + z i 
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P = 



For checking stability, we project any initial perturbation onto the null space of P T , 
where 

ei e 2 Rzx (Sz 1 + 2z 1 )/2 
d e 2 Pz 2 (-3i 2 + 2^ 2 )/2 
ei e 2 Rz 3 (Sz 3 + 2z 3 )/2 
ei e 2 _Rii (— 3ai — ii)/2 
ei e 2 Pi 2 (— 3a 2 — i 2 )/2 
ei e 2 Pi 3 (-3a 3 -i 3 )/2 
The projection matrix is given by 

n = /-p(p T p)- 1 p T . 

From the fact that Zi + z 2 + z 3 = and + i 2 + i 3 = 0, it follows that all columns of P 
are mutually orthogonal except for the 5-th and 6-th columns. Hence, P T P is not a purely 
diagonal matrix. 



Let 



'n-\ 



^n=[Jl(I + hA'(C(U)))). 



,i=0 



We say that an orbit is stable if all eigenvalues of 



lim A„n 

n^oo 



are at most one in magnitude. 



6.1. Stable Orbits 

We computed A n for n — 10 6 . Table 1 shows maximum eigenvalues for those orbits that 
seemed stable from simulation. Table 2 shows maximum eigenvalues for those orbits that 
appeared unstable when simulated. 

Acknowledgements. The author received support from the NSF (CCR-0098040) and 
the ONR (N00014-98-f-0036). 
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Name 


max(Aj(A)) 


max(Aj(AII)) 


Lagrange2 


1.383 


1.362 


FigureEight3 


1.228 


4.220 


Ducati3 


1.105 


3.885 


Hill3_15 


1.444 


2.403 


DoubleDouble5 


12.298 


12.298 


DoubleDoublelO 


1.404 


5.948 


DoubleDouble20 


1.890 


1.890 



Table 1: Apparently stable orbits. 



Name 


max(Aj(A)) 


max(Aj(An)) 


Lagrange3 


81.630 


81.630 


OrthQuasiEllipse4 


18.343 


18.343 


Rosette4 


1.873 


4.449 


Braid4 


727.508 


711.811 


Trefoil4 


41228.515 


41213.852 


FigureEight4 


221.642 


194.095 


FoldedTriLoop4 


74758.355 


74675.092 


PlateSaucer4 


3653.210 


3653.210 


BorderCollie4 


188.235 


188.052 


Trefoil5 


1.913e+8 


1.917e+8 


FigureEight5 


2223.137 


2223.457 



Table 2: Apparently unstable orbits. 
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param N := 3; # number of masses 

param n := 15; # number of terms in Fourier series representation 
param m := 100; # number of terms in numerical approx to integral 

set Bodies := {0. .N-l}; 

set Times := {0..m-l} circular; # "circular" means that next(m-l) = 



param theta {t in Times} := t*2*pi/m; 
param dt := 2*pi/m; 



param aO {i in Bodies} default 
var as {i in Bodies, k in l..n} 
var ac {i in Bodies, k in l..n} 



param bO {i in Bodies} default 0; 
= 0; var bs {i in Bodies, k in l..n} := 0; 
= 0; var be {i in Bodies, k in l..n} := 0; 



var x {i in Bodies, t in Times} 

= a0[i]+sum {k in l..n} ( as [i ,k] *sin(k*theta[t] ) + ac [i ,k] *cos(k*theta[t] ) ); 
var y {i in Bodies, t in Times} 

= b0[i]+sum {k in l..n} ( bs [i ,k] *sin(k*theta[t] ) + be [i ,k] *cos(k*theta[t] ) ); 

var xdot {i in Bodies, t in Times} = (x [i ,next (t)] -x [i ,t] ) /dt ; 
var ydot {i in Bodies, t in Times} = (y [i ,next (t)] -y [i ,t] ) /dt ; 

var K {t in Times} = 0.5*sum {i in Bodies} (xdot[i,t]~2 + ydot [i , t] "2) ; 

var P {t in Times} 

= - sum {i in Bodies, ii in Bodies: ii>i} 

l/sqrt((x[i,t]-x[ii,t])~2 + (y [i ,t] -y [ii ,t] ) "2) ; 



minimize A: sum {t in Times} (K[t] - P[t])*dt; 



let {i in Bodies, k 

let {i in Bodies, k 

let {i in Bodies, k 

let {i in Bodies, k 



in 1 . .n} as [i ,k] := 

in 1 . .n} ac [i ,k] := 

in n. .n} bs [i ,k] : = 

in n. .n} be [i ,k] := 



l*(Uniform01()-0.5) ; 
l*(Uniform01()-0.5) ; 
0.01*(Uniform01O-0.5) ; 
0.01*(Uniform01()-0.5) ; 



solve; 



Fig. 1. — ampl program for finding trajectories that minimize the action functional. 
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param N := 3; # number of masses 

param n := 15; # number of terms in Fourier series representation 

param m := 99; # terms in num approx to integral. must be a multiple of N 

param lagTime := m/N; 

set Bodies := {0. .N-l}; 

set Times := {0..m-l} circular; # "circular" means that next(m-l) = 

param theta {t in Times} := t*2*pi/m; 
param dt := 2*pi/m; 

param aO default 0; param bO default 0; 

var as {k in 1 . . n} : = ; var bs {k in 1 . . n} : = ; 

var ac {k in 1 . . n} : = ; var be {k in 1 . . n} : = ; 

var x {i in Bodies, t in Times} 

= aO+sum {k in l..n} ( as [k] *sin(k*theta[(t+i*lagTime) mod m] ) 

+ ac[k]*cos(k*theta[(t+i*lagTime) mod m] ) ); 
var y {i in Bodies, t in Times} 

= bO+sum {k in l..n} ( bs [k] *sin(k*theta[(t+i*lagTime) mod m] ) 

+ bc[k]*cos(k*theta[(t+i*lagTime) mod m] ) ); 

var xdot {i in Bodies, t in Times} = (x [i ,next (t)] -x [i ,t] ) /dt ; 
var ydot {i in Bodies, t in Times} = (y [i ,next (t)] -y [i ,t] ) /dt ; 

var K {t in Times} = 0.5*sum {i in Bodies} (xdot [i ,t] ~2 + ydot [i , t] "2) ; 

var P {t in Times} 

= - sum {i in Bodies, ii in Bodies: ii>i} 

l/sqrt((x[i,t]-x[ii,t])-2 + (y [i ,t] -y [ii ,t] ) "2) ; 

minimize A: sum {t in Times} (K [t] - P[t])*dt; 



let {k in l..n} as [k] 

let {k in 1 . . n} ac [k] 

let {k in n. .n} bs [k] 

let {k in n. .n} be [k] 

solve ; 



= l*(Uniform01O-0.5) ; 
= l*(Uniform01()-0.5) ; 
= 0.01*(Uniform01()-0.5) ; 
= 0.01*(Uniform01()-0.5) ; 



Fig. 2. — AMPL program for finding choreographies by minimizing the action functional. 
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FigureEight3 Braid4 Trefoil4 




FigureEight4 FoldedTriLoop4 Trefoil5 



FigureEight5 
Fig. 3. — Periodic Orbits — Choreographies. 
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Rosette4 PlateSaucer4 
Fig. 4. — Periodic Orbits-Non-Choreographies. 



BorderCollie4 
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Ducati3_2 Ducati3_0.5 Ducati3_0.1 




Ducati3_alluneq Ducati3_alluneq2 
Fig. 5. — Periodic Orbits — Ducati's with unequal masses. 
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Fig. 6. — Periodic Orbits — Hill-type with equal masses. 



